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Abstract. Let F(x) = X^eN d FvX v be a multivariate power series with complex co- 
efficients that converges in a neighborhood of the origin. Assume F = G/H for some 
functions G and H holomorphic in a neighborhood of the origin. We derive asymptotics 
"^j", for the coefficients F ra as r — > oo with ra <G N d for a in a permissible subset of d-tuples 

of positive reals. More specifically, we give an algorithm for computing arbitrary terms 
of the asymptotic expansion for F ra when the asymptotics are controlled by a transverse 
04 ■ multiple point of the analytic variety H = 0. This improves upon earlier work by R. 

(3JT)[ Pemantle and M. C. Wilson. 

^ ■ We have implemented our algorithm in Sage and apply it to obtain accurate numerical 

results for several rational combinatorial generating functions. 

CO 

O ! 1. Introduction 

In [PW02, PW04] Pemantle and Wilson began a program of analytic combinatorics in 
several variables to derive asymptotic expansions of coefficients of combinatorial gener- 
ating functions. In this article we continue that program by improving upon several of 
their results. 
Let F{x) = 

^ d be a complex power series with complex coefficients 
that converges in a neighborhood Q of the origin. Assume F = G/H for some functions 
G and H holomorphic on Q. For example, F could be a rational function. We derive 
asymptotics for the "ray coefficients" F ra as r — > 00 with ra G N d for a in a permissible 
i/-) ■ subset of d-tuples of positive reals. 

■ In [PW02], Pemantle and Wilson derived the general form of the asymptotic expansion 

of F ra for directions whose asymptotics are controlled by smooth points of the set V = 
{x G f2 : H(x) = 0} of singularities of F, that is, points where V is locally a complex 
manifold. They gave an explicit formula for the leading term but no practical method 
for computing higher order terms. In [RW08], a computational extension of [PW02], we 
devised an algorithm and Maple implementation for computing these higher order terms. 

In [PW04], Pemantle and Wilson generalized [PW02] by deriving the form of the as- 
ymptotic expansion of F ra for directions whose asymptotics are controlled by multiple 
points of V, that is, points where V is locally a finite union of complex manifolds. Again, 
they gave an explicit formula for the leading term but no practical method for computing 
higher order terms. This article is a computational extension of [PW04] analogous to 
[RW08]. Herein we devise an algorithm and Sage [S + f2] implementation for computing 
these higher order terms. 

Why is it important in asymptotic analysis to have algorithms to compute higher order 
terms? There are several reasons. First, the form of the algorithm or formulas itself can 
be insightful. For example, the recasting of smooth point results in terms of Gaussian 
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curvature in [BBBP11] yields a much clearer understanding, independent of coordinates, 
of how local geometry controls the asymptotic scale. Second, computing higher order 
terms often gives one higher numerical accuracy at small values of r than using the leading 
term alone. Third, sometimes computing higher order terms is necessary. For example, 
computing the variance of random variables via a generating function often requires third 
order asymptotics. Fourth, computing higher order terms can be difficult or downright 
infeasible by hand — indeed, this is usually the case with multivariate asymptotics. An 
algorithm allows the end user to pass on the task to a computer. 

Our Contribution. In this article we give an algorithm for computing arbitrary terms 
of the asymptotic expansion for F ra for directions whose asymptotics are controlled by a 
multiple point of V of order n > 1. We do this by first deriving an explicit formula in 
Section 3 for the special case where n < d and the ideal generated by the germ of H in 
the ring of germs of holomorphic functions is radical. This generalizes the formula for 
the smooth point case n = 1 in [RW08, Theorem 3.2] and improves upon the formula in 
[PW04, Theorem 3.5], which gave an explicit formula for only the leading term. We then 
show in Section 5 how to reduce the general multiple point case to the special case. This 
gives a unified method for the computation of higher-order asymptotics that works for any 
value of n and d. Our method of derivation uses Fourier-Laplace integrals as in [PW04], 
but avoids the complications of infinite stationary phase sets. We have implemented 
our algorithm in an open-source Sage file called amgf . sage that is downloadable from 
Raichev's website, and in Section 6 we employ it to work out examples. Section 7 contains 
most of our proofs. 

To the best of our knowledge, the algorithms here and in [RW08] are the first explicit, 
practical, and fairly general methods in the multivariate combinatorics literature for 
computing higher order asymptotic expansions. 

2. Preliminaries 

Throughout this article we make use of basic facts from local analytic geometry, a good 
reference for which is [dJPOO]. 

For brevity we write a power series ^„ gN d a v {x\ — pi) Ul --- (xd — pd) Ud as J2u &v{x — p) v 
and use the multi-index notation v\ = V\\ ■ ■ ■ z/^!, rv = (rvi, . . . , rvd), v + 1 = {y\ + 
1, . . . , Vd + 1), and d u = d^ 1 ■ ■ ■ d^ d , where dj is partial differentiation with respect to 
component j. 

Let 0(Q) denote the C- algebra of holomorphic functions on an open set tt C C d and 
O p the C-algebra of germs of holomorphic functions at p G C d . The latter algebra is a 
local Noetherian factorial ring whose unique maximal ideal is the set {/ G O p : f(p) = 0} 
of non-units. 

We refer often to both <i-tuples and (d — l)-tuples and write a = (a\, . . . , a<z-i) given a 
tuple a = (ai, . . . , ad)- For simplicity we assume d > 2, though our formulas below also 
apply in the case d = 1 of univariate functions, after making the simple changes described 
in [RW08, Remark 3.6]. 

Let Q C C d be a neighborhood of the origin (an open subset of C d containing the origin) 
and F(x) = F v x u G 0(Q). Assume F = G/H for some relatively prime G,H G 0(Q). 
Let V be the set of singularities of F, namely the analytic variety {x G £1 : H(x) = 0} 
determined by H. We will derive asymptotics for the ray coefficients F ra as r — > oo with 
ra G N d for a in a permissible subset of , the set of d-tuples of positive reals. For 
asymptotics of F v when d = 2 and v — > oo along more general paths see [Lla06]. 

To begin we recall several key definitions from [PW02, PW04]. 
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Just as in the univariate case, asymptotics for the coefficients of F are determined by 
the location and type of singularities of F, that is, by the geometry of V. Generally 
the singularities closest to the origin are the most important. We define 'closest to the 
origin' in terms of polydiscs. For p G C d , let D(p) = {x G C d : Vj \xj\ < \pj\} and 
C(p) = {x G C d : Vj \xj\ = \pj\} be the respective polydisc and polycircle centered at the 
origin with polyradius determined by p. 

Definition 2.1. We say that a point p G V is minimal if V D D(p) is contained in the 
boundary of D(p), that is, if there is no point x G V such that for all j, \xj\ < \pj\. We 
say that p G V is strictly minimal if V D D(p) = {p}, and we say that p is finitely 
minimal if V D D(p) is finite. 

Note that V always contains minimal points. To see this, let p G V and define / : 
V fl D{p) — > R by /(a;) = a/o^ + ■ ■ ■ + x\. Since / is a continuous function on a compact 
space, it has a minimum, and that minimum is a minimal point of V. 

The singularities of F with the simplest geometry are the smooth/regular points of V. 
Asymptotics for F ra dependent on smooth points were derived in [PW02, RW08]. Here 
we focus on asymptotics dependent on points with the next simplest geometry, that is, 
multiple points. 

Definition 2.2. Let p G V and consider the unique factorization of the germ of H in 
O p into irreducible germs. Choosing representatives for these germs gives a factorization 
H = H^ 1 ■ ■ ■ H^ n valid in a neighborhood of p. We say that p is a multiple point of 
order n if 

• for all j we have Hj(p) = 0, and 

• every set of at most d vectors from {Vi^i(p), . . . , VH n (p)} is linearly independent. 

We say that p is a convenient multiple point of order n if p is a multiple point of 
order n and there exists an index k such that for all j we have p k dkHj(p) ^ 0. 

In other words, p is a multiple point of V iff V is locally a union of n complex manifolds 
that intersect transversely at p *. In particular, the multiple points of V of order n = 1 
are exactly the smooth points of V, and so multiple points are generalizations of smooth 
points. Notice also that the definition above depends only on information about H in an 
arbitrarily small neighborhood of p and so it is independent of the germ representatives 
chosen. Lastly, to derive an asymptotic expansion of the coefficients F ra we will need 
to consider the singularities of F relevant to the direction a. We call these singularities 
critical points, and they arise when approximating the Fourier-Laplace integrals we use 
to approximate F ra (in Lemmas 4.4 and 4.5). They also have a stratified Morse theoretic 
interpretation which, in the interest of simplicity, we will not pursue here; for more details 
see [PW08, Section 3.1]. 

Definition 2.3. Let a G R^ and let p G V be a convenient multiple point and choose 
an index k such that pkdkHj(p) ^ for all j — 1, . . . , n. Consider the scaled logarithmic 
gradient vectors 



7j(p) = 



p 1 d 1 H j (p) p d d d Hj{p) 
p k d k Hj(p) ' ' " ' p k d k Hj(p) 



) 



*In keeping with [PW08] we are simplifying matters by assuming transversality. For a more general 
definition of 'multiple point' see [PW04]. 
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for j = 1, 



, n. We say that p is critical for a if 



/ a\ a d \ , N 



for some s.,- > 0, that is, if a lies in the conical hull of the ~fj(p), which we call the critical 
cone of p. 

3. The full asymptotic expansion: special case 

Let p G V be a convenient multiple point of order n, and let H = H® 1 ■ ■ ■ H^ n be a 
local factorization of H about p as above. 

Without loss of generality and for concreteness and ease of notation, suppose p d d d Hj{p) ^ 
for all j. Henceforth we breaking symmetry and base our explicit calculations on the 
index d. For instance, when we talk about critical points, we divide by the index-c? terms 
p d d d Hj(p). 

Remark 3.1. For the remainder of this section we assume the special case of a\ = 
. . . = a n = 1 and n < d. 

Now, to state our main results we need to define several auxiliary sets and functions, 
most of which are derived from G and H and arise from the integration tricks we use 
to approximate F ra . The reader should feel free to skim over these definitions on a first 
reading, and move on to the main results starting at Theorem 3.3. 

We parametrize the dth coordinate in terms of the first d — 1 coordinates. Since 
d d Hj(p) ^ for all j, we can apply the Weierstrass preparation theorem to each Hj to 
get 



Hj(w,y) = Uj(w,y) y 



hj(w) 



in a neighborhood of p, where Uj is holomorphic and nonzero at p, hj is holomorphic in 
a neighborhood of p with l/hj(p) = p d , and d d Hj(w, l/hj(w)) ^ 0. Thus 



H(w,y) =U(w,y)\ 



n 

y 



hj (w) 



3=1 

a, n 



- — - / 1 \ 



[since ai 



in a neighborhood of p, where U = U\ ■ ■ ■ U n . We use reciprocals, because they turn out 
to be convenient for proving Lemma 4.5 later on. 

For n > 2 let A = {s e W 1 " 1 : Sj > for all j and Y^Zi s j < !}> the standard 
orthogonal simplex of dimension n — 1. This simplex comes from the residue calculation 
in Lemma 4.2. 

Let W be a neighborhood of p on which the hj are defined. For j — 0, . . . , n — 1 and 
a ERi define the functions h : W x A C, Aj : dom([7) -> C, e : [-1, l]"*" 1 -> C^ 1 , 
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and Aj, h, $ : e'^W PI C(p)) x A -> C by 



= 777 r M 

U(w,y)*-± y 



n-l 



h(w, s) = s-Ji-i H h s n _i/i n _i + (1 - ^2 s i) h r 

3=1 

Aj{w,y) = (-iy-iy-^ (^-) J G(w,y^) 



e(t) = (piexp(iti), . . . ,p d _iexp(it d _ 1 )) 
fc(t,s) = h(e(t),s) 
A j (t,s) = A j (e(t),h(t,s)) 

d-i 

s) = - log(p d /i(t, s))+iV — t m 

m=l 



Note that F(w,y) = G(w, y)J YTj=i{v~ l ~ hj( w )) an d that h, Aj, and $ are all p°° 
functions. The function h comes from the residue calculation in Lemma 4.2, and the 
functions e, h, Aj, and $ come from the exponential change of variables in Lemma 4.3. 

Let J\ og (H,p) denote the n x d logarithmic Jacobian matrix, the jth row of which 
is the logarithmic gradient vector V\ og Hj(p) = (pidiHj(p), . . . ,p d d d Hj{p)). Notice that 
if the convenient multiple point p has all nonzero coordinates, then every subset S C 
{Vi og i?i(p), • • • ,V\ og H n (p)} spans a subspace of C d of dimension \S\. Logarithmic gra- 
dients arise, essentially, from the exponential change of variables used to get a Fourier- 
Laplace integral in Lemma 4.3. 

If a is critical for p, then 



_ . a d s 1 a d s n . 

\p d d d H 1 ( P y- Pd d d H n (p) ] 



for some nonnegative tuple s* with X^=i s j ~ Moreover, if p has all nonzero coordi- 
nates, then the tuple s* is unique since Ji og (p) has rank n < d. Let 9* = (0, ... ,0, si, ... , s*_ 
IR^" 1 x A c W i+n - 2 . ^ 

If the Hessian det <&"(#*) is nonzero, then p is called nondegenerate for a. Critical 
points and nondegeneracy come into play in Lemmas 4.4 and 4.5. 



Remark 3.2. In the smooth point case n = 1 we can simplify the definitions above. In 
that case H = Hf 1 with a\ — 1 (in this section) and we set 



h(w) = hx(w) 
A (w) = y- l G(w,y- 1 ) 

A (t) = A (e(t)) 
hit) = h(e(t)) 



y=h(w) 



d-1 



$(t) = -log(M(t)) + iV^ m 

i Old 
m=l 

e* = t* = o. 

With the setup above we can now get to our main theorem. It is an elaboration of the 
following formula that appeared in [PW04]. 

Theorem 3.3. Let a G and p G V be a strictly minimal convenient multiple point 
with all nonzero coordinates that is critical and nondegenerate for a. Then there exist 
constants b q (a) such that 



F ~ r~ ra 

1 ra ^ 



q>0 



n-d)/2-q 



as r — > oo with ra G N d . Here $"(0*) is the (d + n — 2) x (d + n — 2) Hessian matrix of 
$ evaluated at 



Proof. Proved in [PW04, Theorem 3.9]. 



□ 



We give an explicit formula for all the coefficients bq(pt). Previously, only &o(o0 was 
known. 

Theorem 3.4. Let a G and p G V be a strictly minimal convenient multiple point 
with all nonzero coordinates that is critical and nondegenerate for a. Then 



(*) 



F =r~ ra 

± ra 



N-1 



X 



(27r)("- d )/ 2 (det$"(r))- 1 / 2 ^(ra d ) (n - d)/2 ^ 



0<j <min{n— l,q} 
max{0,g— n}<k<q 
j+k<q 



9=0 






n-j 


C; 1 ) 


n + k — q 



\q-j-k 



as r — > 00 and ra G N . 



6 



Here 



2k "/k+i{~ 



1=0 y ' 
$(#) = $(0) - $ 



U k+l {Aj$){0*) 



2 k + i i\(k + iy. 



0*)$"(0*)(0-0*) 



*\T 



the differential operator TL is given by 

« = - £ 



l<a,fe<d+n-2 



and [?] denotes the Stirling numbers of the first kind. In every term of Lk(Aj,<t) the 

total number of derivatives of Aj and of $" is at most 2k + j. 

Moreover, for each positive integer N the big-oh constant of (*) stays bounded as a 
varies within a compact subset of M.^ of the critical cone of p. 



Proof. In the next section. 



□ 



Remark 3.5. In the smooth point case n — 1, (*) agrees with the formula in [RW08, 
Theorem 3.2]. Moreover, in that case we can allow coordinates of p to be zero as long as 
PkdkH(p) ^ for some k. Also, when n — 1 and d = 2 we can drop the nondegeneracy 
hypothesis ([RW08, Theorem 3.3]). 

Proposition 3.6. Under the hypotheses of Theorem 3.4 we have 

where A is a (d — 1) X (d — 1) matrix, C is an (r — 1) x (d — 1) real matrix, and 

p t diH k (p) pidiH n {p) 



kl 



Pdd d H k (p) p d ddH n (jp) ' 
Notice that we only take derivatives with respect to t in A. 

Proof. Since $ is C°°, its Hessian matrix is symmetric. The formula for A follows by 
definition. To compute the remainder of the Hessian, let s n = J2j< n s i f° r notational 
convenience. For I < d we have 



<9$ 

dti 



(0,8) 



=r-^= h 1 

/i(t,s) «rf 
-ipdPi^Sjdihjjp) + 

-fn(e(t)) + h n (e(t)) 



(0,8) 



. Oil 



h(t,s) 



(0,8) 



0. 



By the implicit function theorem we have d[hj(w) = hj(w) 2 d[Hj(w, 1 / hj(w)) / ddHj(w , 1/hj 
for I < d, j < n, and w G W. So for k,l < d we have 



<9 2 $ 
dskdi. 



(0, s) = -ip d pi(dih k (p) - dih n {p)) (since dsjdsj = -1) 

PidiH k (p) _ pidiH n {p) 
p d d d H k (p) Pdd d H n (p) 



^(0, S )=0, 



as desired. 



Finally each G R, because by [PW08, Proposition 3.12] each 



Pid t H k (p) 
Pdd d H k {p) 



□ 



Theorem 3.7. Under the hypotheses of Theorem 3.4, when n = d there exists e G (0, 1} 
such that 

±G(p) 



p 



+ 0(e r 



det H'(p) U 3 < d Pj 

as r — > oo. Here H'(p) is the n x d Jacobian matrix of if evaluated at p. 

Moreover, the big-oh constant stays bounded as a varies within a compact subset of 
R^ of the critical cone of p. 

Proof. By [PW08, Corollary 3.24] all terms beyond the leading term in the asymptotic 
expansion of p ra F ra are zero and the error term is exponentially decreasing. (This follows 
from a Leray residue argument on the Cauchy integral of F ra .) 



According to (*) the leading term is 



Lq(Aq,$) 
\/dct($"(0*)) 



. First, 



L (A ,<j>)=H'(A )(d*)=A Q (p,-) 

Pd 



(-ir-y d G(p) = (-i)-v^ n ^ 



-G(p) 



-Gjp) _ 

U(p)p d UjKnPdOdHjip) 



Second, by Proposition 3.6, ydet($"(0*)) = ^(det C) 2 = | detC| since C is a real matrix. 
Now consider the n x d matrix T whose jth row is the scaled logarithmic gradient vector 

™ — ( PUhHjM Pdd d Hj(p) 

h \P) I Pd a d Hj ( P ) ' • • • ' Pd a d ff 3 (p) 



. Then 
/ 7i(p) -7„(p) \ 



det C = det 



7n-i(p) - 7n(p) 

\ 7n(p) / 

(by expanding the latter n x d matrix by minors along its last column) 
detr 

(by similarity via elementary row operations) 

Uj<nPj 



UiKnPdddHjip) 



det H\p). 



This proves the result. 



□ 



4. Proving Theorem 3.4 

To prove Theorem 3.4 we follow an approach similar to that of [PW02, PW04, RW08]. 
However, in contrast to those articles, here we first assume that H has the relatively 
simple local factorization H = Hi ■ ■ ■ H n with n < d and then show in Section 5 how to 
reduce to this case. We take the following steps. 

Step 1: Use Cauchy's integral formula to express p ra F ra as a <i-variate integral over 
a contour C in Q. 

Step 2: Expand the contour C across pd and use Cauchy's residue theorem to ex- 
press the innermost integral residue. 
Step 3: Rewrite the residue as an n-variate integral over the simplex A. 
Step 4: Rewrite the resulting integral as a Fourier-Laplace integral. 
Step 5: Approximate the integral asymptotically. 

Starting at step 1, we use Cauchy's integral formula to write 

v rap _ ^P(0) 1 f G(w)d W 

V 1 ra y i ,\ y 



(ra!) * (2vri) d J c w ra + 1 H( W y 

where C is a contour in Q. We then follow steps 2-5 by applying the following lemmas, 
the proofs of which have been swept away to Section 7 to clarify the logical flow of the 
main argument. 

Lemma 4.1 (for step 2). Let a G R+ and p G V be a strictly minimal convenient multi- 
ple point with nonzero coordinates. There exists e G (0, 1) and a polydisc neighborhood 
D of p such that 



f a F ra = f a {2^if- d J ~^dw + O {e r 



as r — > oo with ra G N d , where X = DC\C{p) and R r {w) is the sum over j of the residues 
of y (->■ y~ rad ~ 1 F(u;, y) at hj(w). 

Proof. Proved in [PW02, proof of Lemma 4.1]. □ 



Lemma 4.2 (for step 3). In the previous lemma for n > 2 we have 

ds, 

y=h(w,s) 



d\ n—1 



Rr(w)= j [j-] (-l)-7r(«>,V) 



where f r (w, y) = —y rad 1 G{w, y x ) and ds is the standard volume form ds\ A • • • A ds n _i. 
For the smooth case n = 1 we have R r (w) = f r (w, h(w)). 

Proof. See Section 7. □ 

For j = 0, ... ,n - 1 define Pj : N -> N by P^r) = ( n ~ l )(ra d - l)"" 1 ^'. The falling 

factorial powers in Pj are defined by a- = a(a — 1) ■ • • (a — k + 1) and a- = 1 for a G R 
and fcGE So the degree of Pj in r is n — 1 — j . 

Lemma 4.3 (for step 4). For n > 2, 

re— 1 „ „ 

= (2 7 r) 1 ^^P i (r) / / A,-(t, s) exp(-ra d $(t, s))ds dt + 0(e r ), 



i=o 



9 



as r — ?• oo with ra G N d , where X = e 1 (X). For n — 1, 

p ra F rQ = (27r) 1 - d / 2 J (t)exp(-ra d $(t)) c/t + 0(e r ), 
as r — 7- oo with ra G N d . 

Proof. See Section 7. □ 

The next lemma on Fourier-Laplace integrals provides our key approximation. The 
function spaces mentioned are complex valued. A stationary and nondegenerate point of 
a function g is a point 8* such that Vg(8*) = and det g"(8*) ^ 0, respectively. 

Lemma 4.4 (for step 5). Let £ C M. m be open, iV a positive integer, and g = N+ \m/2~\ . 
If A G C 29 (£) with compact support in £, $ G C 3 « +1 (£), K$ > 0, 5K$(#*) = 0, $ has a 
unique stationary point 6* G supp A, and 0* is nondegenerate, then 

J A(6) exp(-u$(6))d6 = exp(-cj$(r))(det ( ^^ ))~ 1/2 ^"^(A $ )+° ( 
as a; — )• oo. 

Here is the function defined in Theorem 3.4 with m = d + n — 2. Moreover, the 
big-oh constant is bounded when the partial derivatives of $ up to order 3q + 1 and the 
partial derivatives of A up to order 2q all stay bounded in supremum norm over S . 

Proof. Proved in [H6r83, Theorem 7.7.5]. □ 
The final lemma ensures that the hypotheses of Lemma 4.4 are satisfied in our setting. 

Lemma 4.5 (for step 5). Let a G M.+ and p be a strictly minimal convenient multiple 
point that is critical and nondegenerate for a. Then on X x A, we have 3?$ > with 
equality only at points of the form (0, s) (and only at zero for n — 1), and $ has a unique 
stationary point at 9*. 

Proof. See Section 7. □ 

We can now prove Theorem 3.4. 

Proof of Theorem 3.4- By Lemmas 4.1 and 4.3 there exists e G (0, 1) and an open bounded 
neighbourhood X of such that 

n-l 

p ra F ra = (27r) 1 - d ]T Pj(r)I jtr + O (e r ) 

3=0 

as r — ?> oo with ra G N d , where Ij >rn = J £ Aj{6) exp(—rad&(0))d9 and £ = X x A , where 
A° is the interior of A. 

Choose k G C°°(£) with compact support in £ (a bump function) such that k = 1 on 
a neighbourhood Y of 8*. Then 

I j r = I K{8)A j {8)exp{-ra d ${8))d8+ I "(1 - K(8))A j (8)exp(-ra d ^(8))d8. 
Je Je 

The second integral decreases exponentially as r — > oo since 3?$ is strictly positive on 
the compact set £ \ Y by Lemma 4.5. By Lemma 4.5 again and the nondegeneracy 
hypothesis, we we may apply Lemma 4.4 to the first integral. Noting that Lk(nAj, $) = 
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Lk(Aj, $) because the derivatives are evaluated at 9* and k = 1 in a neighborhood of 6* 
we get 



I j>r = exp(-n d $(0*))(det 



ra d §'' 



2tt 



AT-l 
k=0 



N-l 



= (27r)( d+n - 2 )/ 2 (det &'(6*)y l l 2 ^ L fc (A,-, ^(ra,)-^"- 2 '/ 2 -' + O {{ra d )^ d+n - 2 ">/ 2 - N ) 

k=0 

as r — > oo with ra G N d . 

Notice that for j — 0, . . . , n — 1 each Ij T has error 0((rarf)~^ +n_2 ^ 2_Ar ) and each Pj(r) 
has degree r — j — 1 in n. Thus the error in the asymptotic expansion for p ra F ra will be 
a sum of terms of the form 0{{ra d ) {n - d)/2 - N - j ) which is 0{{ra d f n - d)l2 - N ). So 

N-1 

p ra F ra = (2vr)( ri - d )/ 2 (det ^'(fl*))" 1 / 2 ^ ^(a)(r« d ) (n - d)/2 " 9 + O ((ro^)^/ 2 ^) 

9=0 

n-1 N-1 

= (27r)^- d )/ 2 (det &\6*))- 1 ' 2 £ ^ P,(r)L fc (I„ S)(ra d )-( d+n - 2 >/ 2 - fc + O ((m,,)-^-^ . 

j=0 fc=0 



Let us expand fj(r) and collect like powers to find the coefficients b q (a) 
The falling factorial powers satisfy (a — 1)— = (a — 1) ... (a — 1 — k) = 
related to regular powers and Stirling numbers of the first kind via 



-a^^- and are 



m r n 

m 



1=0 



■l) m - l a l - 



see [GKP94, (6.13)] for instance. Thus 



n — 1 



£ 

z=o 



n — 3 
I 



{-l) n -i-\ra d ) 1 



and so 

AT-l 



^ b q (a)(ra d ) 

9=0 



n-1 AT-l n— j 
j=0 fc=0 i=0 





"n - j" 


CT) 





'-^"^'-'(TOd)-^/ 2 -^. 



The coefficient 6 9 (a) is found by imposing the constraint (n—d)/2 — q = —(d+n)/2 — k + l. 
Thus I = n + k — q, and we can eliminate the /-sum to arrive at formula (*). 

Lastly, regarding uniformity, we may assume that the Aj and $ are defined and hence 
C°° on a neighborhood of the closure of £, so that their derivatives up to any given order 
all stay bounded in supremum norm over £. Now suppose a varies within a compact 
subset K C of the critical cone of p. Since J\ og (H,p) has rank n < d it is a bijective 
linear transformation from W 1 to its image in R d and therefore a bicontinuous function. 
Thus its inverse maps K to a compact set i^' of #*s in S. Choose the neighborhood Y in 
the argument above to contain K' so that one bump function k works for all 6*. Since 
the derivatives of the nAj and $ up to any given order all stay bounded in supremum 
norm over S and since only $ and $' depend on a but continuously, we conclude by 
Lemma 4.4 that for any given N, the big-oh constant in (*) remains bounded as a varies 
within K. □ 

n 



5. The full asymptotic expansion: general case 



Again let p e V be a strictly minimal convenient multiple point of order n with all 
coordinates nonzero and let H = H® 1 ■ ■ ■ H® n be a local factorization of H. We deal now 
with the case of arbitrary aj and n. 

In step 2 of the previous section the Cauchy integral can be manipulated to reduce to 
the special case a\ = . . . = a n = 1 and n < d. More specifically, we amend our plan by 
inserting these three steps after step 2: 

(2a) If r > d, then decompose F as a sum of fractions whose denominators are of type 
Y\j £ j H- 3 where J is a size d subset of {1, . . . ,n} and each bj is an integer with 
bj < aj. So each denominator in the sum has only d irreducible factors of H. 

(2b) If some irreducible factor of H is repeated, then treat each resulting integral 
as the integral of a holomorphic form, and rewrite each integral as the sum of 
integrals whose denominators are of type w ra+1 Y\j € j Hj where J is a size at most 
d subset of {1, . . . , n}. So each holomorphic form has a denominator with at most 
d unrepeated irreducible factors of H . 
(6) Add up all the asymptotic expansions. 

The following two lemmas prove that these additional steps are possible. 

Lemma 5.1 (for step 2a). Let p be a multiple point of H :— H® 1 ■ ■ ■ H% H , where n > d, 
each Hj is holomorphic in a neighborhood U of p, and the germ of each Hj is prime. 
Then for any function G holomorphic on U, there exists a neighborhood of p in which we 
have the partial fraction decomposition 

G_ _ \ ^ Gj 
H~ ^\\ #V 

where each Gj is holomorphic (and possibly zero), J ranges over all subsets of {1, . . . , n} 
of size d, and for each J we have J2jej bj = Y?2=i a «- 

Proof. Since p is a multiple point of H, the gradients at p of any d of the Hj are linearly 
independent. Thus the germs of any d of the Hj generate the maximal ideal in O v by 
[Rui93, Corollary 5.4]. In particular, the germ of H\ is in the ideal of the germs of 
H 2 , . . . , H n , and so in a neighborhood of p we have 

n 

i=2 

for some holomorphic functions gj. Therefore, in that neighborhood we have 

G_ G Y,'j=2 9jHj 
H H1 1+1 H% 2 ■■■HZ- 

° 2 + 9i + + 



H ai+l H a 2 -l H a 3 . . . ffa n Hf 1+1 Hf 1 H? • • • H*> 
Gjn 

where Gj = Ggj. Notice that in the denominator of each resulting summand, the sum of 
the degrees of all the Hj remains Y17=i ai - 

Recursively repeating this procedure on each summand (always singling out Hi, say) 
yields the desired result in finitely many steps. □ 
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Remark 5.2. When the Hj are polynomials from a computable polynomial ring, such 
as Q[x], the procedure in the proof above is computable. Alternatively a partial fraction 
expansion can be computed according to the algorithm in [Lei78], which is not applicable 
to the analytic case. 

Lemma 5.3 (for step 2b). Let p be a multiple point of H := H^ 1 ■ ■ ■ H n n , where n < d, 
each Hj is holomorphic in a neighborhood U of p, and the germ of each Hj is prime, and 
let Vj := {x G U : Hj(x) = 0}. Then for any function G holomorphic on U, there exists 
a neighborhood U' of p such that the holomorphic form 

H(x) 

is de Rham cohomologous in U' \ (Vi U • • • U V n ) to a holomorphic form 

Gj(x)dx 

where each Gj is holomorphic (and possibly zero) on U' and J ranges over all subsets 
of {1, ... ,n}. In particular, the integrals of the two forms above over a polycircle in 
U' \ (Vi U ■ • • U V„) are equal. 

Proof. Proved in [AY83, Theorem 17.6]. □ 

Remark 5.4. When the Hj are polynomials from a computable polynomial ring, such 
as Q[x], the procedure in the proof above is computable. 

When applying Lemma 5.3 in step 2b to our local integrals of residues, G(x) will be 
of the form G(x)/x ra+l where G(x) does not contain r and for each J we will have 
^2j eJ bj = Xir=i a «- Thus upon inspection of the constructive proof of Lemma 5.3, the 
cohomologous form will have r-degree at most YH=i( a i ~~ -0> w here the powers of n arise 
from the derivatives of G(x)/x ra+1 . 

In particular, if n > d and the other assumptions of Theorem 3.7 hold, then we can 
combine Lemmas 5.1 and 5.3 and Theorem 3.7 to conclude that the leading term of the 
asymptotic expansion of p ra F ra is a polynomial of degree at most Y^i=i a * ~ as ^ s a ^ so 
shown in [PW04, Theorem 3.6]. 

Remark 5.5 (for step 6). When computing the asymptotics for G/H in a direction a 
by summing up the asymptotic contributions from the terms of the form Gj/Y\j eJ Hj 
where J has size at most d, the only terms that will contribute to the expansion (modulo 
an exponentially decreasing error term) are the ones whose critical cone (the conical hull 
of {jj(p) '■ j G J}) contains a and whose numerator does not vanish at p [PcmOO, Section 
5]. In the case where all such contributing terms have numerators that vanish at p, a 
finer analysis is required to determine the correct asymptotics of G/H which we do not 
provide here (but will be included in Pemantle and Wilson's forthcoming book on analytic 
combinatorics in several variables). 

Remark 5.6. In case p is finitely minimal, for each point x of V D C(jp) we simply find 
an open set around x and apply the general procedure above. After that we sum the 
resulting asymptotic expansions over the finitely many x. 
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6. Examples 



Let us apply the formulas and procedures of Sections 3 and 5 to a few combinatorial 
examples, that is, to functions with all nonnegative Maclaurin coefficients. We will use 
our Sage package amgf . sage. 

We focus on combinatorial examples F(x), because for any a G M.+ there is a minimal 
point in Vfl that determines the asymptotics for F ra ([PW08, Theorem 3.16]). 

Since there is no known computable procedure to factor an arbitrary polynomial H in 
the analytic local ring of germs of holomorphic functions about p, we choose examples 
where if is a polynomial whose local factorization in the algebraic local ring about p 
equals its factorization in the analytic local ring about p, that is, H is a polynomial 
whose irreducible factors in C[x] are all smooth at p. 

Example 6.1 {n < d, no repeated factors). Consider the trivariate rational function 

F(x,y,z) 



{l-x(l + y))(l-zx 2 (l + 2y)) 

in a neighborhood Q of the origin; cf [PW08, Example 4.10]. Its coefficients F u are all 
nonnegative, and its denominator H(x,y,z) factors over C[x, y,z] into irreducible terms 



x(l + y) and H 2 (x, y, z) 



zx 



'1 + 2y), both of which are globally 



Hi(x,y,z) = 
smooth. 

The set of non-smooth/singular points of V = {(x,y,z) G f2 : H(x,y,z) = 0} is 
V = {(x, y, z) G Q : H(x, y, z) = VH(x, y, z) = 0} = {(l/(a + 1), a, (a + 1) 2 /(2a + 1)) : 
a G C \ { — 1}}, which consists entirely of convenient multiple points of order n = 2. They 
are not convenient in coordinate d = 3, but are in coordinate d — 1 =2, which we use here 
for our calculations. A simple check shows that the points (l/(a + 1), a, (a + 1) 2 / (2a + 1)) 
for a > are strictly minimal. 

The critical cone for each such point is the conical hull of the vectors 71 = (1, a/ (a + 
1),0) and 72 = (l,a/(2a + l),l/2). 

For instance, p = (1/2, 1,4/3) controls asymptotics for all a in the conical hull of the 
vectors 71 (p) = (2,1,0) and ^(p) = (3,1,3/2). For instance, a = (8,3,3) is in this 
critical cone, and applying Theorem 3.4 we get 

1231 329047 



108 r 



r -l/2 _ r -3/2 + °^^ ±l r -5/2 + 0( r -7/2) 



_V21^ 8232V21T' ' 58084992 v^vr 

as r —7- 00. 

Calling the one-term, two-term, and three-term truncations of this asymptotic formula 
Si (r), 5 2 (r), and 5 3 (r), respectively and comparing them with the actual values of F ra 
for small r, we get the following table. 



r 


1 


2 


4 


8 


108" r F TO 




0.3518518519 


0.2548010974 


0.1823964231 


0.1297748629 


108- r Si(r) 




0.3693487820 


0.2611690282 


0.1846743909 


0.1305845142 


108- r S 2 (r) 




0.3509381749 


0.2546598957 


0.1823730650 


0.1297708726 


108- r S 3 (r) 




0.3516356189 


0.2547831876 


0.1823948602 


0.1297747255 


108" r 5i(r) rel 


err 


-0.04972811712 


-0.02499177148 


-0.01248910347 


-0.006238891584 


108- r S 2 (r) rel 


err 


0.002596766210 


0.0005541644108 


0.0001280622701 


0.00003074786527 


108- r S 3 (r) rel 


err 


0.0006145569473 


0.00007028933620 


0.000008568698736 


0.000001058756657 



Table 1. Successive approximations to p 
a = (8,3,3). 



x F rn with relative errors for 



14 



Example 6.2 (n < d, no repeated factors). Consider the trivariate rational function 



v ' y ' ' (A-2x-y-z){A-x-2y-z) 

in a neighborhood Q of the origin; cf [PW04, Example 3.10]. Its coefficients F v are all 
nonnegative, and its denominator H(x, y, z) factors over C[x, y, z] into irreducible terms 
Hi(x,y,z) = 4 — 2x — y — z and H 2 (x, y,z) = 4 — x — 2y — z, both of which are globally 
smooth. 

The set of non-smooth points of V = {(x, y, z) G fl : H(x, y, z) = 0} is V' = {(x, y, z) G 
VL : H (x, y, z) = VH(x, y,z) = 0} = {(1 — a, 1 — a, 1 + 3a : a G C}, which contains a 
line segment {(1 — a, 1 — a, 1 + 3a : —1/3 < a < 1} of convenient multiple points of order 
n = 2. The convenient multiple point p = (1,1,1) is strictly minimal and its critical cone 
is the conical hull of the vectors ji(p) = (2, 1, 1) and ^(p) = (1, 2, 1). 

For instance, a = (3, 3, 2) is in the critical cone and applying Theorem 3.4 we get 

F - =jm ( 4r " /2 - rr 3 ' 2 + W'"' 2 ) + 

as r — )• oo. 

Calling the one-term, two-term, and three-term truncations of this asymptotic formula 
Si(r), S 2 (r), and S 3 (r), respectively and comparing them with the actual values of F ra 
for small r, we get the following table. 



r 


1 


2 


4 


8 


16 


± rOL 






0.7849731445 


0.7005249476 


0.5847732654 


0.4485547669 


0.3237528587 


Si(r) 






1.302940032 


0.9213177319 


0.6514700159 


0.4606588663 


0.3257350080 


S 2 (r) 






1.189837598 


0.8813299831 


0.6373322117 


0.4556603976 


0.3239677825 


S 3 (r) 






1.202663729 


0.8835973440 


0.6377330283 


0.4557312524 


0.3239803079 


Si(r) 


rel 


err 


-0.6598530037 


-0.3151819005 


-0.1140557451 


-0.02698466340 


-0.006122414820 


S 2 (r) 


rel 


err 


-0.5157685415 


-0.2580993527 


-0.08987918808 


-0.01584116640 


-0.0006638514355 


Ss(r) 


rel 


err 


-0.5321081198 


-0.2613360125 


-0.09056461026 


-0.01599912872 


-0.0007025396085 



Table 2. Successive approximations to p ra F ra with relative errors for 
a = (3,3,2). 



Notice that in this case the three-term approximation to F ra is not an improvement 
over the two-term approximation for n < 16. The question, which we do not discuss 
here, of how many terms of a divergent asymptotic series expansion to use for a given 
argument to obtain the best approximation/least error is called the question of 'optimal 
truncation' or 'optimal approximation'. See [PK01], for instance, for more details. 

Example 6.3 (n < d, repeated factors). Consider the trivariate rational function 



v ' y ' ' (4-2x-y-z) 2 (4:-x-2y-z) 

in a neighborhood Q of the origin. Its coefficients F v are all nonnegative, and its denom- 
inator H (x, y,z) = (4 — 2x — y — z) 2 (A — x — 2y — z) is shown factored over C[x, y, z\. 
Since H contains repeated factors, we first reduce 

F(x, y, z) dx A dy A dz 

rgCtin+l y02n+l gCtsn+l ' 
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the differential form of the Cauchy integral of F, to a de Rham cohomologous form with 
no repeated factors, namely 

[16(2a 3 y — ct2z)n + 16(2?/ — z)]/ (yz) dx A dy A dz 

(4 - 2x - y - Z){A - X - 2y - Z ) x <xin+l ya2 n+l z a 3 n+l ' 

which determines the asymptotics of F ra . The constructive proofs of Lemma 5.1 (in the 
case of polynomials) and Lemma 5.3 to find such a cohomologous form are implemented 
in amgf . sage. 

The singular variety V of this new form is the same as in the previous example and 
so the singularity analysis is the same. The convenient multiple point p = (1,1,1) is 
strictly minimal and its critical cone is the conical hull of the vectors ji(p) = (2, 1, 1) and 
72 (p) = (1,2,1). 

Taking a = (3, 3, 2) again, for instance, and applying Theorem 3.4 we get 

F = — ( Ar 1 ' 2 + — r" 1 ' 2 - 1967 r -V* \ + 0(r^' 2 ) 
ra VS^V 72 41472 ) +U{ h 

as r — > oo. 

It is a coincidence that the leading coefficient above is the same as the leading coefficient 
in the previous example without repeated factors. Using the denominator (4 — 2x — y — 
z) 3 (4 — x — 2y — z) instead, for instance, gives a different leading coefficient. 

Calling the one-term, two-term, and three-term truncations of this asymptotic formula 
Si(r), 5*2 (r), and Ss(r), respectively and comparing them with the actual values of F ra 
for small r, we get the following table. 



r 


1 


2 


4 


8 


16 


TP 

± ra 






0.9812164307 


1.576181132 


2.485286378 


3.700576827 


5.260983954 


Si(r) 






1.302940032 


1.842635464 


2.605880063 


3.685270927 


5.211760127 


S2(r) 






1.515572607 


1.992989400 


2.712196350 


3.760447895 


5.264918270 


S 3 (r) 






1.500123128 


1.987527184 


2.710265167 


3.759765118 


5.264676873 


Si(r) 


rel 


err 


-0.3278824031 


-0.1690505784 


-0.04852305395 


0.004136084917 


0.009356391776 


S2(r) 


rel 


err 


-0.5445854345 


-0.2644418586 


-0.09130133815 


-0.01617884746 


-0.0007478289298 


S 3 (r) 


rel 


err 


-0.5288402039 


-0.2609763838 


-0.09052429168 


-0.01599434190 


-0.0007019445473 



Table 3. Successive approximations to p ra F ra with relative errors for 
a = (3,3,2). 



Notice that in this case the two-term or three-term approximation to F ra is not an 
improvement over the one-term approximation until somewhere between r = 8 and r = 
16. 

Example 6.4 (n > d with no repeated factors). Consider the bivariate function 



v { i_2x-y){l-x-2y) 

in a neighborhood Q of the origin; cf [PW08, Example 4.12]. 

Its coefficients F v are all nonnegative, and its denominator H(x,y) factors over C[x, y] 
into irreducible terms Hi(x, y) = 1 — 2x — y and H 2 (x,y) = 1 — x — 2y, both of which are 
globally smooth. 

The set of non-smooth points of V = {(x,y) G f2 : H(x,y) = 0} is V = {(x,y) G Q : 
H(x, y) = VH(x, y) = 0}, which consists of the convenient multiple point p = (1/3, 1/3) 
of order r = 2. The point p is strictly minimal and its critical cone is the conical hull of 
the vectors 71 (p) = (2, 1) and 72 (p) = (1/2, 1). 
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By Theorem 3.7, there exists e G (0, 1) such that for any a in this critical cone we get 

F « = p_ra ( u ( fL ,i + 0(e r )) = 3^+ Q2 >"(3 + 0(e r )), 

VP1P2I det#'(p)| J 

as r —7- 00. 

Taking a = (4, 3), say, letting S(r) be the asymptotic expansion above, and comparing 
it to the actual values of F ra for small r, we get the following table. 



r 


1 


2 


4 


8 


16 


2187 r F ra 
2187- r S(r) 
2187" r 5'(r) rel err 


1.960219479 
3 

-0.5304408677 


2.298399383 
3 

-0.3052561804 


2.587511051 
3 

-0.1594153382 


2.809909562 
3 

-0.06765002002 


2.950100341 
3 

-0.01691456340 



Table 4. Successive approximations to p ra F ra with relative errors for 
a = (4,3). 



Example 6.5 (n > d with repeated factors). Consider the bivariate function 

= {l-2x-yf{l-x-2y) 2 ' 

which is a variation of the function of the previous example. 

Since the denominator of F contains repeated factors, we first reduce 

F(x, y) dx A dy 

the differential form of the Cauchy integral of F to a de Rham cohomologous form with 
no repeated factors, which amgf . sage computes. 

Reusing the analysis of the previous example and applying Theorem 3.7, there exists 
e G (0, 1) such that for any a in conical hull of the vectors 71 (p) = (2, 1) and 72 (p) = 
(1/2,1) we get 

F ra = 3 (ai+Q2)n (-3(2a? - 5«ia 2 + 2a 2 2 )r 2 - 3(«i + a 2 )n - 9 + 0(e r )), 

as r — > 00 and for some e G (0, 1). 

Taking a = (4, 3), say, letting S(r) be the asymptotic expansion above, and comparing 
it to the actual values of F ra for small r, we get the following table. 



r 


1 


2 


4 


8 


16 


2187 r Fra 


30.72702332 


111.9315678 


442.7813138 


1799.879232 


7367.545085 


2187- r S(r) 





69 


387 


1743 


7335 


2187- r S"(r) rel err 


1.000000000 


0.3835519207 


0.1259793763 


0.03160169385 


0.004417358124 



Table 5. Successive approximations to p ra F ra with relative errors for 
a = (4,3). 



17 



7. Remaining Proofs 
Proof of Lemma 4-2. Let f r (w,y) = —y rad ~ l G(w , y^ 1 ) . Then for r > 2, 



i?» = V lim j/"™"^ - ^H- 1 )^^, </) 

™ 41 i. fl ■ I 111 1 1 



. 1 y->hj(w)- 



Jfcr- -iT^M-V - ftj ^^ ni.i(y 1 -M»)) 

f r (w,hj(w)) 

c/(7i y dcr 2 • • • y ( 7^ J / r (w, (1 - (7i)/ii + (cri - CT 2 )/i2 H 

(cr r _ 2 - cr n _i)/i n _i + a n -ih n ) da n ^i 

(by [DL93, Chapter 4, Section 7, equations (7.7) and (7.12)]) 

/. / Q \ n—1 
J \j^J (-l) n - 1 f r (w,S 1 h 1 + --- + S n . 1 h n _ 1 + (l-J2s j )h n )ds 

(by the change of variables (si, . . . , s n -i) = (1 — ci> Ci — 02, . . . , o>_2 — <r n _i)), 



as desired. 

Notice that the (— l)"" 1 cancels with the (— l)™^ 1 in the definition of f r . 

For n = 1, we have i? r (w) = lim^^)-! y~ rad ~ l (y - h (w)~ 1 )F(w, y) = f r (w,h(w)). 

□ 



Proof of Lemma 4-3. First, for n > 2, 



<9y 



(-ir-7(^,y) 

n-l 



i=o 

n-l 

i=o 
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Thus 



x w r&+1 



=p ra (2vri) 1 - d / -L- f 
Jxw ra+1 J A 

(by Lemma 4.2) 



d\ n— 1 



ds dw 

y=h(w,s) 



9y 



n-1 



j=0 

n-1 



= (2 7 ri) 1 -^^W / / ^*>, h(w, s))(p d h(*>, s)) rad ds *™ 

j=Q JxJaW \[ m =l W rn 

n— 1 „ „ d— 1 

=(27r) 1 - d VP.W / / Uexp(-ia m ntJA 3 (t,s)(p d h(t,s)Y ad dsdt 

j=0 Jx J A m=1 

(via the change of variables w = e(t)) 

n—l „ „ 

=(27r) 1_d y] p j( r ) L / Aj(t, s) exp(-ra d $(t, 8))ds dt, 

which with Lemma 4.1 proves the stated formula for p ra F ra . 

The formula for the case n = 1 follows similarly. □ 

Proof of Lemma 4-5. First $(0, s) = and 

n 

m(t,s) = - log \p d h(t,s)\ > -log^s,|Mi(e(t))| > 

j=i 

for t 7^ 0, because the sum is convex and \hj{w)~ 1 \ > \p d \ for w ^ V since p is strictly 
minimal. 

Now by the calculation in the proof of Proposition 3.6, for all / < d we have 

9 1 $(n = -iE44§M+^ = ". 

ipdddHjip) a d 

where the last equality holds since p is critical for a. Also di$>(6*) = for d < I < n+d—2 
since $(0, s) is constant. Thus V$(6 1 *) = 0. Now det $"(0*) 7^ 0, since p is nondegenerate 
for a. So there is a neighborhood of 9* in which 6** is the only zero of Thus, shrinking 
X x A if needed, #* is the unique stationary point of $. □ 
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